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Abstract 

We determine the topological susceptibility xt m two-flavor QCD using the lattice simulations 
at a fixed topological sector. The topological charge density is unambiguously defined on the 
lattice using the overlap-Dirac operator which possesses exact chiral symmetry. Simulations are 
performed on a 16 3 x 32 lattice at lattice spacing ~ 0.12 fm at six sea quark masses m q ranging 
in m s /6-m s with m s the physical strange quark mass. The Xt is extracted from the constant 
behavior of the time-correlation of flavor-singlet pseudo-scalar meson two-point function at large 
distances, which arises from the finite size effect due to the fixed topology. In the small m q 
regime, our result of xt is proportional to m q as expected from chiral effective theory. Using the 
formula xt = m g E/iV/ by Leutwyler-Smilga, we obtain the chiral condensate in Nf = 2 QCD 
as S MS (2 GeV) = [252(5)(10)MeV] 3 , in good agreement with our previous result obtained in the 
e-regime. 

PACS numbers: 11.15.Ha,11.30.Rd,12.38.Gc 
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The vacuum of Quantum Chromodynamics (QCD) has a non-trivial topological structure. 
The cluster property and the gauge invariance require that the ground state must be the 
9 vacuum, a superposition of gauge configurations in different topological sectors. The 
topological structure is also essential for the U(l) problem, i.e. the flavor singlet pseudo- 
scalar meson knows the presence of non-trivial topological charge in the QCD vacuum. 
Therefore, the topological excitations, such as the instantons, played a central role in the 
past attempts to understand the QCD vacuum since the early years of QCD. 

The topological charge fluctuations are characterized by the topological susceptibility Xt, 
defined as Xt = (Q 2 ) through the topological charge Q in a given volume Q. Since the 
topological excitations do not occur in the perturbation theory, theoretical calculation of 
Xt from the QCD Lagrangian necessarily involves non-perturbative methods, such as the 
numerical simulation of lattice QCD. 

Lattice calculation of Xt has become reasonably precise in the quenched approximation 
, which is however not realistic as the quark loop effects are discarded. On the other 
hand, despite much effort, more realistic calculation without quenching remains a challenging 
problem. For instance, in full QCD, xt is expected to linearly depend on the quark mass and 
to vanish in the chiral limit. Even such a basic property has not been clearly reproduced in 
the previous lattice calculations (for a recent study, see, e.g., J3J). 

The main difficulties are the following, (i) Definition of the topological charge den- 
sity using the gauge links (discretization of (l/32ii 2 )e l j iU p a tr(F l j il/ F pa )) causes bad ultraviolet 
divergences, while the cooling methods devised to tame the short distance fluctuations intro- 
duce sizable systematic uncertainties, (ii) Unquenched simulations with Wilson/staggered 
fermion do not respect correct chiral or flavor symmetry at finite lattice spacing, and the 
definition of the topological charge through the Atiyah-Singer index theorem is ambiguous, 
(iii) With the molecular-dynamics-type algorithms, which are based on a continuous evolu- 
tion of the gauge links, the system is trapped in a fixed topological sector as the continuum 
limit is approached. Therefore, a proper sampling of different topological sectors cannot be 
achieved. Approaching the chiral limit (or the physical up and down quark masses) also 
makes the tunneling of topological charge a rare event, because of the suppression of the 
fermion determinant for large topological charges. 

During the last decade, (i) and (ii) have been solved by the realization of exact chiral 
symmetry on the lattice, with which the topological charge is uniquely defined at any finite 
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lattice spacing by counting the number of fermionic zero-modes. (For recent studies re- 
specting the exact chiral symmetry, we refer .) However, (iii) remains insurmountable, 
since the correct sampling of topology becomes increasingly more difficult towards realistic 
simulation with lighter quarks and finer lattices. A plausible solution is to perform QCD 
simulations in a fixed topological sector and to extract Xt from local topological fluctua- 
tions. In j(| (see also {7]), a general formula to transcribe any observable measured at a 
fixed topological charge to its value in the 9 vacuum is derived. As an application, a new 
method to calculate Xt is also proposed. 

In this paper, we use this method to precisely calculate Xt in two-flavor lattice QCD 
with exact chiral symmetry. The results are compared with the prediction from chiral 
perturbation theory (%PT) jg] 

Xt = ^ + 0(m 2 q ), (1) 

for Nf flavors of sea quarks with mass m q . The chiral condensate £ can be determined 
independently, e.g. from the low-lying eigenvalues of the Dirac operator [h]]. Therefore, 
the comparison provides a critical test of the lattice approach to study the QCD vacuum in 
the chiral regime. 

A two-point function of the topological charge density p(x) calculated in a finite volume 
VL at a given topological charge Q behaves as [6] 

^J p{xM0))Q= ^- Xt -^j + O^) (2) 

where C4 = — ((Q 4 ) — 3(Q 2 ) 2 )/Q. The expectation value (• • -)q denotes an average in a given 
topological sector Q. The correlation does not vanish even for large separations, because of 
the violation of the clustering property at fixed topological charge. We emphasize that the 
derivation of (j2J) relies only on modest assumptions such as (Q 2 ) ^> 1 and Q <C (Q 2 ), which 
are the conditions to apply the saddle point expansion in the Fourier transform from a fixed 
9 to a fixed Q. Except for these conditions, the formula is model independent. 

We consider, in particular, two spatial sub- volumes at t\ and t2, for which the correlator 
is defined as 

C(ti - h) = (g(ti)Q(fe)) = E {/>(xi)pM) , (3) 

Xl,&2 

where the summations run over the spatial sites x± and £2 at t± and £2, respectively. Its 
plateau at large \ti — t 2 | can be used to extract Xt, provided that \ci\ 2x^(1 
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FIG. 1: A schematic diagram for the time-correlation function of the flavor singlet operator P°(x). 
Each solid line denotes the valence quark propagator. 



In order to preserve the exact chiral symmetry, which is essential for the definition of the 
topological charge, we employ the overlap-Dirac operator [12I, Q 



Tfl \ f Tfl 

D(m q ) = (m + -^-j + ( m - -±\ 7 5 sgn [H w (-m Q )} (4) 

with mass m q . The kernel operator H\y{— mo) is the conventional Wilson-Dirac operator 
with a large negative mass term — m . 

In place of the topological charge density p(x) (and p(0)) in ([2]), we use m q P°(x) (and 
m q P°(0)), that were shown to give the same asymptotic constant as (J2J) p] (the origi- 
nal suggestion is in [ll|), where P°(x) is the flavor singlet pseudo-scalar density P°(x) = 
^E%$ f {x)i 5 [l-aD(0)/(2mo)]il) f (x). The correlator C^t) = (P°(x)P°(0)) contains 
a connected and a disconnected diagram as shown in Fig. [TJ If we pick the disconnected piece 
and identify a "topological charge density" , it can be written as pi(x) = m g tr[7 5 (D c +m (? )~J.], 
where D c is a chirally-symmetric (y^Dr + -D c 7s = 0) nonlocal operator, relating to -D(O) by 
D c = [1 — aD(0)/(2m )]^ 1 D(0) [14|. Integrated over the entire lattice volume, p\(x) reduces 
to the number of fermionic zero-modes, and thus has the necessary property for the topolog- 
ical charge density. This implies that the correlator (p 1 (a;)p 1 (0)) has the same asymptotic 
constant as ([2]). However, the correlator (m q P°(x)m q P°(0)) approaches the constant with 
the rate governed by the r( mass, eT m v'\ x \ which is much faster than ^- m ^\ x \ appearing in 
(pi(x)pi(O)). 

Simulations are carried out for two-flavor (Nf = 2) QCD on a 16 3 x 32 lattice at a lattice 
spacing ~ 0.12 fm. For the gluon part, the Iwasaki action is used at (3 = 2.30 together with 



unphysical Wilson fermions and associated twisted-mass ghosts 



151 ] . The unphysical degrees 



of freedom generate a factor det[Hyy(— rao)/{Hy^{— mo) + /i 2 )] in the partition function (we 
take to = 1.6 and \i = 0.2) that suppresses the near-zero eigenvalue of Hw(— TO)) and thus 
makes the numerical operation with the overlap operator (J3J) substantially faster. Further- 
more, since the exact zero eigenvalue is forbidden, the global topological change is preserved 
during the molecular dynamics evolution of the gauge field. Our main runs are performed 
at Q = 0, while Q = —2 and —4 configurations are also generated at one sea quark mass 
in order to check the consistency as described below. Ergordicity within a given global 
topological charge is satisfied if the configuration space of that topological sector forms a 
connected manifold. This is indeed the case in the continuum SU(3) gauge theory on a 
four- dimensional torus, and therefore is probably also true at small lattice spacing adopted 
in this work. 

We use the Hybrid Monte Carlo algorithm 



with the mass preconditioning 



13| . The 



fermion masses for the preconditioner were chosen to be 0.4 for heavier sea quark masses and 
0.2 for the two lightest ones (see later). We exploit the rational approximation a la Zolotarev 
for the sign function in (00) after projecting out low-lying eigenmodes of Hy/{— too). With 
the number of poles in the rational function to be 8-10, the accuracy of O(10~^ 7_8 - ) ) is 
achieved for the sign function. The simulations have been done in two phases for each sea 
quark mass. In the first phase the nested conjugate gradient (CG) is used to invert the 
overlap operator (jlj) (see [l8,Q for details). On the other hand, in the second phase we use 
the five-dimensional implementation of the overlap solver without the low-mode projection. 
The target accuracy of O(10^ 7-8 )) is maintained by adding an additional Metropolis step 
calculated with the nested CG [201 ] . 

For the sea quark mass m q we take six values: 0.015, 0.025, 0.035, 0.050, 0.070, and 
0.100 that cover the mass range m s /6-m s with m s the physical strange quark mass. After 
discarding 500 trajectories for thermalization, we accumulate 10,000 trajectories in total for 
each sea quark mass. In the calculation of Xt, we take one configuration every 20 trajec- 
tories, thus we have 500 configurations for each m q . For each configuration, 50 conjugate 
pairs of lowest-lying eigenmodes of the overlap-Dirac operator D(0) are calculated using the 
implicitly restarted Lanczos algorithm and stored for the later use. In these calculations, 
the better accuracy of O(10~ 12 ) is enforced for the sign function by increasing the number 
of poles in the rational approximation. 

For the connected diagram (see Fig. [1]), the pion correlator is computed using the con- 
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FIG. 2: Time-correlation function of the flavor singlet rj' (circles) and its connected (triangle 
down) and disconnected (triangle up) contributions. Data at m q = 0.025 are shown. 

jugate gradient algorithm with a low-mode preconditioning. Low-modes are also used for 



averaging over source points 



2 ll ]. which significantly improves the statistical signal. For the 



disconnected diagram, the quark propagator is represented by the eigenmode decomposi- 
tion and approximated by the 50 conjugate pairs of the low-lying eigenmodes. The quark 
propagator is then obtained for any source point without extra computational cost, and the 
disconnected loops can be calculated with an average over the source point. The trunca- 
tion is motivated by the expectation that the long distance correlation is dominated by the 
low-lying fermion modes; its validity has to be checked numerically (see below). 

In Fig. [21 we plot Cy(i) together with those of connected and disconnected parts for 
m q = 0.025. The curve is the fit to the function A + B{er Mt + e~ M ( T ~*)) with data for 
C n >(t) in the range t G [4, 28]. The horizontal line is the fitted constant A, and the curvature 
of the data points represents the contamination by the flavor-singlet state rj', that rapidly 
decays due to its heavy mass. Assuming |c 4 | <C 2x^fi, we obtain a A \t = 3.40(27) x 10~ 5 
at m q = 0.025. The error is estimated using the jackknife method with bin size of 20 
configurations, with which the statistical error saturates. 

Since the disconnected diagram is computed with only 50 pairs of low-lying eigenmodes, 
we have to check whether they suffice to saturate Cy (t) . For the time range [4, 28] used for 
fitting, as the number of eigenmodes is increased from 10 to 30, the change of correlator 
[SCrfl/Crf is ~ 3%, while from 30 to 50, it is only ~ 0.3%, which is less than 8% of the 
statistical error. Thus C v > is well saturated with 50 eigenmodes. This also holds for all six 
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FIG. 3: Topological susceptibility xt r o versus sea quark mass m q ro. 
sea quark masses. 

In Fig. El we plot the topological susceptibility Xt r o as a function of the sea quark mass 
rnqTo, where ro is the Sommer scale extracted from the static quark potential. We also 
present our data in Table fl] together with the lattice spacing a determined from the static 



quark potential with the input r = 0.49 fm [19], and the pseudoscalar mass m^ro [22] . 
The statistical precision is good enough to find a clear dependence on the sea quark mass. 
For the smallest three quark masses, 0.015, 0.025, and 0.035, the data are well fitted by a 
linear function F + Gm with the intercept F = 0.2(3) x 10~ 3 and the slope G = 0.087(4). 
Evidently, the intercept is consistent with zero, in agreement with the %PT expectation (CQ). 
Equating the slope to r^Ti/Nf, we obtain r^E = 0.174(8). 

In order to convert E to that in the MS scheme, we calculate the renormalization factor 
Z^f{2 GeV) using the non-perturbative renormalization technique through the RI/MOM 



scheme 



23fl. Our result is Z™ s (2GeV) = 0.742(12) [22j With a" 1 = 1670(20)(20) MeV 
determined with r = 0.49 fm Q], the value of X is transcribed to S M5 (2 GeV) = 
252(5) (10) MeV] 3 , which is in good agreement with our previous result [251(7) (11) MeV]" 



9|, lid] obtained in the e-regime from the low-lying eigenvalues. The errors represent a 
combined statistical error (or 1 and Z^f) and the systematic error estimated from the higher 
order effects {e.g., C4 term), respectively. Since the calculation is done at a single lattice 
spacing, the discretization error cannot be quantified reliably, but we do not expect much 
larger error because our lattice action is free from 0(a) discretization effects. 

In principle, Xt m (El) is universal for any topological sector. We check the universality of 
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TABLE I: The topological susceptibility x* r o versus the quark mass m q ro obtained in this work. 
The lattice spacing a for each quark mass is determined from the static quark potential with the 
input r = 0.49 fm [19]. 



m q a 


a[fm] 


m n r 


m q r 




x t n 


0.015 


0.1194(15) 


0.7096(102) 


0.0616(8) 


5.59(47) x 10~ 3 


2.582(88) 


0.025 


0.1206(18) 


0.8935(137) 


0.1016(15) 


9.26(1.29) x 10~ 3 


4.454(356) 


0.035 


0.1215(15) 


1.053(13) 


0.1412(17) 


0.0123(22) 


6.078(794) 


0.050 


0.1236(14) 


1.238(14) 


0.1982(22) 


0.0187(34) 


9.900(1.354) 


0.070 


0.1251(13) 


1.456(15) 


0.2742(28) 


0.0263(33) 


14.65(1.21) 


0.100 


0.1272(12) 


1.734(17) 


0.3852(36) 


0.0199(33) 


11.82(1.55) 



Xt as follows. At sea quark mass m q = 0.050, we generate 250 configurations with Q = — 2 
and —4 respectively in addition to the main run at Q — 0. Then we extract Xt from 
the time-correlation function of if, similar to the Q = case. Our results for a 4 Xt are: 
{7.4(1.3), 6.4(2.1), 5.9(1.8)} x 10~ 5 for Q = {0, -2, -4} respectively. Evidently, \t extracted 
from different topological sectors are consistent with each other within the statistical error. 

Finally, we discuss on the neglected terms in ([2]). The leading correction comes from the 
| I 2xt & term. With the formulas derived in 6j, we can obtain an estimate of the upper 
bound of |c 4 |/(2x^f2) by measuring the two-point correlator (pi(xi)pi(x 2 )) and the four- 
point correlator {pi(xi)pi(x2)pi(x^)pi(x^)) . Our preliminary result is |c4 1 / (2%^ SI) < 0.1, for 
all six sea quark masses. Details of this calculation will be presented elsewhere. We note 
that this upper bound is also consistent with the ratio \c^\/xt — 0.3 obtained in the quenched 
approximation 24], [25]. 

In this letter, we have determined the topological susceptibility Xt i n two-flavor QCD 
from a lattice calculation of two-point correlators at a fixed global topological charge. The 
expected sea quark mass dependence of Xt from is clearly observed with the good 

statistical precision we achieved, in contrast to the previous unquenched lattice calcula- 
tions. Our result indicates that the topologically non-trivial excitations (e.g., instanton and 
anti-instanton pairs) are in fact locally active in the QCD vacuum, even when the global 
topological charge is kept fixed. The information of these topological excitations is carried 
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by low-lying fermion eigenmodes if the exact chiral symmetry is respected on the lattice. 
This work demonstrates that Monte Carlo simulation of lattice QCD with fixed topology is 
a viable approach, to be pursued when the topology change hardly occurs near the contin- 
uum limit even with chirally non-symmetric lattice fermions. The artifacts due to the fixed 
topology in a finite volume can be removed to obtain the physics results in the 9 vacuum, 
provided that Xt has been determined in the first place (], Q] as has been done in this work. 
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